DEG analysis after resequencing of the 3 last batches, excluding Nitrogen, LFC=1, including T0
Thomas Dobrenel and Nicolas Delhomme
2019-12-28
- 1 Setup
- 2 Raw data
- 3 Expression data
- 4 Data normalisation
- 5 Multivariate analysis
- 6 Differential expression based on all conditions compared to the T0
- 7 Comparisons of GOI lists (Venn diagrams)
- 8 Export list of DEGs
- 9 GO analysis
- 10 Session Info
1 Setup
- Set the working dir
setwd("/mnt/picea/projects/arabidopsis/jhanson/arabidopsis-nutrition-TOR/T0")- Load libraries
suppressPackageStartupMessages({
library(data.table)
library(DESeq2)
library(gplots)
library(here)
library(hyperSpec)
library(limma)
library(LSD)
library(magrittr)
library(matrixStats)
library(parallel)
library(pander)
library(plotly)
library(RColorBrewer)
library(scatterplot3d)
library(tidyverse)
library(tximport)
library(VennDiagram)
library(vsn)
})- Source some helper functions
suppressPackageStartupMessages({
source("~/Git/UPSCb/UPSCb-common/src/R/featureSelection.R")
source("~/Git/UPSCb/UPSCb-common/src/R/gopher.R")
source("~/Git/UPSCb/UPSCb-common/src/R/plot.multidensity.R")
source("~/Git/UPSCb/UPSCb-common/src/R/volcanoPlot.R")
})- Cutoff Different from Schurch et al., RNA, 2016 (LFC at 1 instead of 0,5)
lfc <- 1
FDR <- 0.01- Create palettes
pal <- c(brewer.pal(8,"Dark2"),1)
pal2 <- brewer.pal(9,"Paired") #require package RColorBrewer
cols <- rainbow(17)
hpal <- colorRampPalette(c("blue","white","red"))(100)- Register the default plot margin
mar <- par("mar")2 Raw data
2.1 Loading of sample information
- Read the sample information
samples <- read.csv("~/Git/UPSCb/projects/arabidopsis-nutrition-TOR/doc/samples3.csv")- Remove unnecessary samples
samples %<>% filter(!grepl("P11554_1",SciLifeID)) %>%
filter(! SciLifeID %in% c("P13406_101",
"P14066_128",
"P14066_133",
"P13406_102",
"P14066_131")) %>%
filter(! Nutrition == "PKS") %>%
mutate(Nutrition,Nutrition=relevel(Nutrition,"NPS")) %>%
mutate(AZD,AZD=relevel(AZD,"DMSO"))
samples <- samples[order(samples$Timepoint, samples$Nutrition, samples$AZD),]
samples %<>% mutate(Timepoint,Timepoint=factor(paste0("T",Timepoint)))2.1.1 Load the data
- Call the data
filenames <- list.files("../Salmon",
recursive = TRUE,
pattern = "quant.sf",
full.names = TRUE)- Name the data
names(filenames) <- sub("_S.*","",sapply(strsplit(filenames, "/"), .subset, 3))- Match data <=> sample list
filenames <- filenames[match(samples$SciLifeID,names(filenames))]
filenames <-filenames[!is.na(names(filenames))]
samples <- samples[match(names(filenames),samples$SciLifeID),]- Annotate the samples
samples$Conditions <- factor(paste(samples$Timepoint,samples$Nutrition,samples$AZD,sep="_"))
samples$Batch <- factor(substr(samples$SciLifeID,1,8))3 Expression data
Read the expression at the transcript level
tx <- suppressMessages(tximport(files = filenames,
type = "salmon",
txOut = TRUE))summarise to genes
tx2gene <- data.frame(TXID=rownames(tx$counts),
GENEID=sub("\\.[0-9]+","",rownames(tx$counts)))
gx <- summarizeToGene(tx,tx2gene=tx2gene)## summarizing abundance## summarizing counts## summarizing lengthkg <- round(gx$counts) Sanity check
stopifnot(all(colnames(kg) == samples$SciLifeID))3.1 Raw data export
#dir.create(file.path("analysis_Tom","salmon"),showWarnings=FALSE,recursive=TRUE)
#write.table(kg,file="/mnt/picea/projects/arabidopsis/jhanson/arabidopsis-nutrition-TOR/analysis_Tom/nutrition-unormalised-gene-expression_data.csv")
#save(kg, samples, file = "/mnt/picea/projects/arabidopsis/jhanson/arabidopsis-nutrition-TOR/analysis_Tom/counts.rda")3.2 Preliminary validations
3.2.1 Check for the genes that are never expressed
sel <- rowSums(kg) == 0
sprintf("%s%% (%s) of %s genes are not expressed",
round(sum(sel) * 100/ nrow(kg),digits=1),
sum(sel),
nrow(kg))## [1] "17.9% (5793) of 32310 genes are not expressed"4 Data normalisation
For visualization, the data is submitted to a variance stabilization transformation using DESeq2. The dispersion is estimated independently of the sample type
dds <- DESeqDataSetFromMatrix(
countData = kg,
colData = samples,
design = ~ Conditions)## converting counts to integer modedds <- estimateSizeFactors(dds)4.1 Perform a Variance Stabilizing Transformation for plotting
vst <- varianceStabilizingTransformation(dds,blind=FALSE)
vsd <- assay(vst)
vsd <- vsd - min(vsd)Write out
#write.csv(vst,"/mnt/picea/projects/arabidopsis/jhanson/arabidopsis-nutrition-TOR/analysis_Tom/library-size-normalized_variance-stabilized_data_nutrition.csv")5 Multivariate analysis
5.1 PCA
Principal Component Analysis on the normalized data * Establishment of the PCA
conditions1 <- factor(paste(samples$AZD,samples$Nutrition,sep="_"))
pc <- prcomp(t(vsd))
percent <- round(summary(pc)$importance[2,]*100);percent## PC1 PC2 PC3 PC4 PC5 PC6 PC7 PC8 PC9 PC10 PC11 PC12 PC13 PC14 PC15 PC16
## 43 24 13 6 4 2 1 1 1 1 0 0 0 0 0 0
## PC17 PC18 PC19 PC20 PC21 PC22 PC23 PC24 PC25 PC26 PC27 PC28 PC29 PC30 PC31 PC32
## 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
## PC33 PC34 PC35 PC36 PC37 PC38 PC39 PC40 PC41 PC42 PC43 PC44
## 0 0 0 0 0 0 0 0 0 0 0 0- Graphical representation of PC1 x PC2
conds <- droplevels(conditions1)
plot(pc$x[,1],
pc$x[,2],
xlab=paste("Comp. 1 (",percent[1],"%)",sep=""),
ylab=paste("Comp. 2 (",percent[2],"%)",sep=""),
col=pal[as.integer(conds)],
pch=c(15,16,17)[as.factor(samples$Timepoint)],
main="All timepoints")
legend("bottomright",pch=19,
col=pal[1:nlevels(conds)],
legend=levels(conds))
legend("topleft",pch=c(15,16,17),
col="black",
legend=c("T0","T6","T24"))
- Graphical representation of PC2 x PC3
plot(pc$x[,2],
pc$x[,3],
xlab=paste("Comp. 2 (",percent[2],"%)",sep=""),
ylab=paste("Comp. 3 (",percent[3],"%)",sep=""),
col=pal[as.integer(conds)],
pch=c(15,16,17)[as.factor(samples$Timepoint)],
main="All timepoints")
legend("topleft",pch=19,
col=pal[1:nlevels(conds)],
legend=levels(conds))
legend("bottomleft",pch=c(15,16,17),
col="black",
legend=c("T0","T6","T24"))
- Graphical representation of PC1 x PC3
plot(pc$x[,1],
pc$x[,3],
xlab=paste("Comp. 1 (",percent[1],"%)",sep=""),
ylab=paste("Comp. 3 (",percent[3],"%)",sep=""),
col=pal[as.integer(conds)],
pch=c(15,16,17)[as.factor(samples$Timepoint)],
main="All timepoints")
legend("bottomleft",pch=19,
col=pal[1:nlevels(conds)],
legend=levels(conds))
legend("bottomright",pch=c(15,16,17),
col="black",
legend=c("T0","T6","T24"))
6 Differential expression based on all conditions compared to the T0
6.1 Preparation of the design
sel <- samples$Timepoint %in% c("T0","T6","T24")
suppressMessages(dds <- DESeqDataSetFromMatrix(
countData = kg[,sel],
colData = samples[sel,],
design = ~ Conditions))6.2 Differential expression analysis
dds <- DESeq(dds)## estimating size factors## estimating dispersions## gene-wise dispersion estimates## mean-dispersion relationship## final dispersion estimates## fitting model and testing## -- replacing outliers and refitting for 2 genes
## -- DESeq argument 'minReplicatesForReplace' = 7
## -- original counts are preserved in counts(dds)## estimating dispersions## fitting model and testing6.3 Variance Stabilising Transformation
6.3.1 Perform a Variance Stabilizing Transformation for plotting
vst <- varianceStabilizingTransformation(dds,blind=FALSE)
vsd <- assay(vst)
vsd <- vsd - min(vsd)6.4 Contrasts to T0
The contrast by default is the first one (not Intercept)
resultsNames(dds)## [1] "Intercept" "Conditions_T24_NP_AZD_vs_T0_T0_0"
## [3] "Conditions_T24_NP_DMSO_vs_T0_T0_0" "Conditions_T24_NPS_AZD_vs_T0_T0_0"
## [5] "Conditions_T24_NPS_DMSO_vs_T0_T0_0" "Conditions_T24_NS_AZD_vs_T0_T0_0"
## [7] "Conditions_T24_NS_DMSO_vs_T0_T0_0" "Conditions_T6_NP_AZD_vs_T0_T0_0"
## [9] "Conditions_T6_NP_DMSO_vs_T0_T0_0" "Conditions_T6_NPS_AZD_vs_T0_T0_0"
## [11] "Conditions_T6_NPS_DMSO_vs_T0_T0_0" "Conditions_T6_NS_AZD_vs_T0_T0_0"
## [13] "Conditions_T6_NS_DMSO_vs_T0_T0_0"6.4.1 Different effects at T6
6.4.1.1 T6_NPS_DMSO
res <- results(dds,name = "Conditions_T6_NPS_DMSO_vs_T0_T0_0")
data <- data.frame(rownames(res), res$log2FoldChange, res$padj)
write.table(data, file = "Conditions_T6_NPS_DMSO_vs_T0_T0_0.csv", sep = ";", row.names = FALSE, dec=",")
cutoffs <- abs(res$log2FoldChange) >= lfc & ! is.na(res$padj) & res$padj <= FDR
cutoff1 <- res$log2FoldChange >= lfc & ! is.na(res$padj) & res$padj <= FDR
cutoff2 <- res$log2FoldChange <= -lfc & ! is.na(res$padj) & res$padj <= FDR
T6NPSDMSO <- rownames(res[cutoffs,])
T6NPSDMSO_Low <- rownames(res[cutoff2,])
T6NPSDMSO_High <- rownames(res[cutoff1,])
message(sprintf("There are %s genes that are differentially expressed",sum(cutoffs)))## There are 5181 genes that are differentially expressedmessage(sprintf("There are %s genes that are induced",sum(cutoff1)))## There are 3158 genes that are inducedmessage(sprintf("There are %s genes that are repressed",sum(cutoff2)))## There are 2023 genes that are repressed6.4.1.2 MA plot
DESeq2::plotMA(res)
6.4.1.3 Volcano plot
volcanoPlot(res, lfc=1)
6.4.1.4 Heatmap
hm <- heatmap.2(t(scale(t(vsd[cutoffs,]))),col=hpal,
Colv=FALSE,dendrogram = "row",trace = "none",labRow = FALSE,
distfun = pearson.dist, labCol = paste(samples$Nutrition,samples$AZD,sep="-")[sel],
hclustfun = function(X){hclust(X,method="ward.D2")},margins = c(6,5))
6.4.1.5 T6_NPS_AZD
res <- results(dds,name = "Conditions_T6_NPS_AZD_vs_T0_T0_0")
data <- data.frame(rownames(res), res$log2FoldChange, res$padj)
write.table(data, file = "Conditions_T6_NPS_AZD_vs_T0_T0_0.csv", sep = ";", row.names = FALSE, dec=",")
cutoffs <- abs(res$log2FoldChange) >= lfc & ! is.na(res$padj) & res$padj <= FDR
cutoff1 <- res$log2FoldChange >= lfc & ! is.na(res$padj) & res$padj <= FDR
cutoff2 <- res$log2FoldChange <= -lfc & ! is.na(res$padj) & res$padj <= FDR
T6NPSAZD <- rownames(res[cutoffs,])
T6NPSAZD_Low <- rownames(res[cutoff2,])
T6NPSAZD_High <- rownames(res[cutoff1,])
message(sprintf("There are %s genes that are differentially expressed",sum(cutoffs)))## There are 4557 genes that are differentially expressedmessage(sprintf("There are %s genes that are induced",sum(cutoff1)))## There are 2619 genes that are inducedmessage(sprintf("There are %s genes that are repressed",sum(cutoff2)))## There are 1938 genes that are repressed6.4.1.6 MA plot
DESeq2::plotMA(res)
6.4.1.7 Volcano plot
volcanoPlot(res, lfc=1)
6.4.1.8 Heatmap
hm <- heatmap.2(t(scale(t(vsd[cutoffs,]))),col=hpal,
Colv=FALSE,dendrogram = "row",trace = "none",labRow = FALSE,
distfun = pearson.dist, labCol = paste(samples$Nutrition,samples$AZD,sep="-")[sel],
hclustfun = function(X){hclust(X,method="ward.D2")},margins = c(6,5))
6.4.1.9 T6_NS_DMSO
res <- results(dds,name = "Conditions_T6_NS_DMSO_vs_T0_T0_0")
data <- data.frame(rownames(res), res$log2FoldChange, res$padj)
write.table(data, file = "Conditions_T6_NS_DMSO_vs_T0_T0_0.csv", sep = ";", row.names = FALSE, dec=",")
cutoffs <- abs(res$log2FoldChange) >= lfc & ! is.na(res$padj) & res$padj <= FDR
cutoff1 <- res$log2FoldChange >= lfc & ! is.na(res$padj) & res$padj <= FDR
cutoff2 <- res$log2FoldChange <= -lfc & ! is.na(res$padj) & res$padj <= FDR
T6NSDMSO <- rownames(res[cutoffs,])
T6NSDMSO_Low <- rownames(res[cutoff2,])
T6NSDMSO_High <- rownames(res[cutoff1,])
message(sprintf("There are %s genes that are differentially expressed",sum(cutoffs)))## There are 3402 genes that are differentially expressedmessage(sprintf("There are %s genes that are induced",sum(cutoff1)))## There are 1916 genes that are inducedmessage(sprintf("There are %s genes that are repressed",sum(cutoff2)))## There are 1486 genes that are repressed6.4.1.10 MA plot
DESeq2::plotMA(res)
6.4.1.11 Volcano plot
volcanoPlot(res, lfc=1)
6.4.1.12 Heatmap
hm <- heatmap.2(t(scale(t(vsd[cutoffs,]))),col=hpal,
Colv=FALSE,dendrogram = "row",trace = "none",labRow = FALSE,
distfun = pearson.dist, labCol = paste(samples$Nutrition,samples$AZD,sep="-")[sel],
hclustfun = function(X){hclust(X,method="ward.D2")},margins = c(6,5))
6.4.1.13 T6_NS_AZD
res <- results(dds,name = "Conditions_T6_NS_AZD_vs_T0_T0_0")
data <- data.frame(rownames(res), res$log2FoldChange, res$padj)
write.table(data, file = "Conditions_T6_NS_AZD_vs_T0_T0_0.csv", sep = ";", row.names = FALSE, dec=",")
cutoffs <- abs(res$log2FoldChange) >= lfc & ! is.na(res$padj) & res$padj <= FDR
cutoff1 <- res$log2FoldChange >= lfc & ! is.na(res$padj) & res$padj <= FDR
cutoff2 <- res$log2FoldChange <= -lfc & ! is.na(res$padj) & res$padj <= FDR
T6NSAZD <- rownames(res[cutoffs,])
T6NSAZD_Low <- rownames(res[cutoff2,])
T6NSAZD_High <- rownames(res[cutoff1,])
message(sprintf("There are %s genes that are differentially expressed",sum(cutoffs)))## There are 4804 genes that are differentially expressedmessage(sprintf("There are %s genes that are induced",sum(cutoff1)))## There are 2695 genes that are inducedmessage(sprintf("There are %s genes that are repressed",sum(cutoff2)))## There are 2109 genes that are repressed6.4.1.14 MA plot
DESeq2::plotMA(res)
6.4.1.15 Volcano plot
volcanoPlot(res, lfc=1)
6.4.1.16 Heatmap
hm <- heatmap.2(t(scale(t(vsd[cutoffs,]))),col=hpal,
Colv=FALSE,dendrogram = "row",trace = "none",labRow = FALSE,
distfun = pearson.dist, labCol = paste(samples$Nutrition,samples$AZD,sep="-")[sel],
hclustfun = function(X){hclust(X,method="ward.D2")},margins = c(6,5))
6.4.1.17 T6_NP_DMSO
res <- results(dds,name = "Conditions_T6_NP_DMSO_vs_T0_T0_0")
data <- data.frame(rownames(res), res$log2FoldChange, res$padj)
write.table(data, file = "Conditions_T6_NP_DMSO_vs_T0_T0_0.csv", sep = ";", row.names = FALSE, dec=",")
cutoffs <- abs(res$log2FoldChange) >= lfc & ! is.na(res$padj) & res$padj <= FDR
cutoff1 <- res$log2FoldChange >= lfc & ! is.na(res$padj) & res$padj <= FDR
cutoff2 <- res$log2FoldChange <= -lfc & ! is.na(res$padj) & res$padj <= FDR
T6NPDMSO <- rownames(res[cutoffs,])
T6NPDMSO_Low <- rownames(res[cutoff2,])
T6NPDMSO_High <- rownames(res[cutoff1,])
message(sprintf("There are %s genes that are differentially expressed",sum(cutoffs)))## There are 4907 genes that are differentially expressedmessage(sprintf("There are %s genes that are induced",sum(cutoff1)))## There are 3050 genes that are inducedmessage(sprintf("There are %s genes that are repressed",sum(cutoff2)))## There are 1857 genes that are repressed6.4.1.18 MA plot
DESeq2::plotMA(res)
6.4.1.19 Volcano plot
volcanoPlot(res, lfc=1)
6.4.1.20 Heatmap
hm <- heatmap.2(t(scale(t(vsd[cutoffs,]))),col=hpal,
Colv=FALSE,dendrogram = "row",trace = "none",labRow = FALSE,
distfun = pearson.dist, labCol = paste(samples$Nutrition,samples$AZD,sep="-")[sel],
hclustfun = function(X){hclust(X,method="ward.D2")},margins = c(6,5))
6.4.1.21 T6_NP_AZD
res <- results(dds,name = "Conditions_T6_NP_AZD_vs_T0_T0_0")
data <- data.frame(rownames(res), res$log2FoldChange, res$padj)
write.table(data, file = "Conditions_T6_NP_AZD_vs_T0_T0_0.csv", sep = ";", row.names = FALSE, dec=",")
cutoffs <- abs(res$log2FoldChange) >= lfc & ! is.na(res$padj) & res$padj <= FDR
cutoff1 <- res$log2FoldChange >= lfc & ! is.na(res$padj) & res$padj <= FDR
cutoff2 <- res$log2FoldChange <= -lfc & ! is.na(res$padj) & res$padj <= FDR
T6NPAZD <- rownames(res[cutoffs,])
T6NPAZD_Low <- rownames(res[cutoff2,])
T6NPAZD_High <- rownames(res[cutoff1,])
message(sprintf("There are %s genes that are differentially expressed",sum(cutoffs)))## There are 3835 genes that are differentially expressedmessage(sprintf("There are %s genes that are induced",sum(cutoff1)))## There are 2311 genes that are inducedmessage(sprintf("There are %s genes that are repressed",sum(cutoff2)))## There are 1524 genes that are repressed6.4.1.22 MA plot
DESeq2::plotMA(res)
6.4.1.23 Volcano plot
volcanoPlot(res, lfc=1)
6.4.1.24 Heatmap
hm <- heatmap.2(t(scale(t(vsd[cutoffs,]))),col=hpal,
Colv=FALSE,dendrogram = "row",trace = "none",labRow = FALSE,
distfun = pearson.dist, labCol = paste(samples$Nutrition,samples$AZD,sep="-")[sel],
hclustfun = function(X){hclust(X,method="ward.D2")},margins = c(6,5))
6.4.2 Different effects at T24
6.4.2.1 T24_NPS_DMSO
res <- results(dds,name = "Conditions_T24_NPS_DMSO_vs_T0_T0_0")
data <- data.frame(rownames(res), res$log2FoldChange, res$padj)
write.table(data, file = "Conditions_T24_NPS_DMSO_vs_T0_T0_0.csv", sep = ";", row.names = FALSE, dec=",")
cutoffs <- abs(res$log2FoldChange) >= lfc & ! is.na(res$padj) & res$padj <= FDR
cutoff1 <- res$log2FoldChange >= lfc & ! is.na(res$padj) & res$padj <= FDR
cutoff2 <- res$log2FoldChange <= -lfc & ! is.na(res$padj) & res$padj <= FDR
T24NPSDMSO <- rownames(res[cutoffs,])
T24NPSDMSO_Low <- rownames(res[cutoff2,])
T24NPSDMSO_High <- rownames(res[cutoff1,])
message(sprintf("There are %s genes that are differentially expressed",sum(cutoffs)))## There are 6256 genes that are differentially expressedmessage(sprintf("There are %s genes that are induced",sum(cutoff1)))## There are 3062 genes that are inducedmessage(sprintf("There are %s genes that are repressed",sum(cutoff2)))## There are 3194 genes that are repressed6.4.2.2 MA plot
DESeq2::plotMA(res)
6.4.2.3 Volcano plot
volcanoPlot(res, lfc=1)
6.4.2.4 Heatmap
hm <- heatmap.2(t(scale(t(vsd[cutoffs,]))),col=hpal,
Colv=FALSE,dendrogram = "row",trace = "none",labRow = FALSE,
distfun = pearson.dist, labCol = paste(samples$Nutrition,samples$AZD,sep="-")[sel],
hclustfun = function(X){hclust(X,method="ward.D2")},margins = c(6,5))
6.4.2.5 T24_NPS_AZD
res <- results(dds,name = "Conditions_T24_NPS_AZD_vs_T0_T0_0")
data <- data.frame(rownames(res), res$log2FoldChange, res$padj)
write.table(data, file = "Conditions_T24_NPS_AZD_vs_T0_T0_0.csv", sep = ";", row.names = FALSE, dec=",")
cutoffs <- abs(res$log2FoldChange) >= lfc & ! is.na(res$padj) & res$padj <= FDR
cutoff1 <- res$log2FoldChange >= lfc & ! is.na(res$padj) & res$padj <= FDR
cutoff2 <- res$log2FoldChange <= -lfc & ! is.na(res$padj) & res$padj <= FDR
T24NPSAZD <- rownames(res[cutoffs,])
T24NPSAZD_Low <- rownames(res[cutoff2,])
T24NPSAZD_High <- rownames(res[cutoff1,])
message(sprintf("There are %s genes that are differentially expressed",sum(cutoffs)))## There are 4367 genes that are differentially expressedmessage(sprintf("There are %s genes that are induced",sum(cutoff1)))## There are 2093 genes that are inducedmessage(sprintf("There are %s genes that are repressed",sum(cutoff2)))## There are 2274 genes that are repressed6.4.2.6 MA plot
DESeq2::plotMA(res)
6.4.2.7 Volcano plot
volcanoPlot(res, lfc=1)
6.4.2.8 Heatmap
hm <- heatmap.2(t(scale(t(vsd[cutoffs,]))),col=hpal,
Colv=FALSE,dendrogram = "row",trace = "none",labRow = FALSE,
distfun = pearson.dist, labCol = paste(samples$Nutrition,samples$AZD,sep="-")[sel],
hclustfun = function(X){hclust(X,method="ward.D2")},margins = c(6,5))
6.4.2.9 T24_NS_DMSO
res <- results(dds,name = "Conditions_T24_NS_DMSO_vs_T0_T0_0")
data <- data.frame(rownames(res), res$log2FoldChange, res$padj)
write.table(data, file = "Conditions_T24_NS_DMSO_vs_T0_T0_0.csv", sep = ";", row.names = FALSE, dec=",")
cutoffs <- abs(res$log2FoldChange) >= lfc & ! is.na(res$padj) & res$padj <= FDR
cutoff1 <- res$log2FoldChange >= lfc & ! is.na(res$padj) & res$padj <= FDR
cutoff2 <- res$log2FoldChange <= -lfc & ! is.na(res$padj) & res$padj <= FDR
T24NSDMSO <- rownames(res[cutoffs,])
T24NSDMSO_Low <- rownames(res[cutoff2,])
T24NSDMSO_High <- rownames(res[cutoff1,])
message(sprintf("There are %s genes that are differentially expressed",sum(cutoffs)))## There are 1794 genes that are differentially expressedmessage(sprintf("There are %s genes that are induced",sum(cutoff1)))## There are 680 genes that are inducedmessage(sprintf("There are %s genes that are repressed",sum(cutoff2)))## There are 1114 genes that are repressed6.4.2.10 MA plot
DESeq2::plotMA(res)
6.4.2.11 Volcano plot
volcanoPlot(res, lfc=1)
6.4.2.12 Heatmap
hm <- heatmap.2(t(scale(t(vsd[cutoffs,]))),col=hpal,
Colv=FALSE,dendrogram = "row",trace = "none",labRow = FALSE,
distfun = pearson.dist, labCol = paste(samples$Nutrition,samples$AZD,sep="-")[sel],
hclustfun = function(X){hclust(X,method="ward.D2")},margins = c(6,5))
6.4.2.13 T24_NS_AZD
res <- results(dds,name = "Conditions_T24_NS_AZD_vs_T0_T0_0")
data <- data.frame(rownames(res), res$log2FoldChange, res$padj)
write.table(data, file = "Conditions_T24_NS_AZD_vs_T0_T0_0.csv", sep = ";", row.names = FALSE, dec=",")
cutoffs <- abs(res$log2FoldChange) >= lfc & ! is.na(res$padj) & res$padj <= FDR
cutoff1 <- res$log2FoldChange >= lfc & ! is.na(res$padj) & res$padj <= FDR
cutoff2 <- res$log2FoldChange <= -lfc & ! is.na(res$padj) & res$padj <= FDR
T24NSAZD <- rownames(res[cutoffs,])
T24NSAZD_Low <- rownames(res[cutoff2,])
T24NSAZD_High <- rownames(res[cutoff1,])
message(sprintf("There are %s genes that are differentially expressed",sum(cutoffs)))## There are 3437 genes that are differentially expressedmessage(sprintf("There are %s genes that are induced",sum(cutoff1)))## There are 1466 genes that are inducedmessage(sprintf("There are %s genes that are repressed",sum(cutoff2)))## There are 1971 genes that are repressed6.4.2.14 MA plot
DESeq2::plotMA(res)
6.4.2.15 Volcano plot
volcanoPlot(res, lfc=1)
6.4.2.16 Heatmap
hm <- heatmap.2(t(scale(t(vsd[cutoffs,]))),col=hpal,
Colv=FALSE,dendrogram = "row",trace = "none",labRow = FALSE,
distfun = pearson.dist, labCol = paste(samples$Nutrition,samples$AZD,sep="-")[sel],
hclustfun = function(X){hclust(X,method="ward.D2")},margins = c(6,5))
6.4.2.17 T24_NP_DMSO
res <- results(dds,name = "Conditions_T24_NP_DMSO_vs_T0_T0_0")
data <- data.frame(rownames(res), res$log2FoldChange, res$padj)
write.table(data, file = "Conditions_T24_NP_DMSO_vs_T0_T0_0.csv", sep = ";", row.names = FALSE, dec=",")
cutoffs <- abs(res$log2FoldChange) >= lfc & ! is.na(res$padj) & res$padj <= FDR
cutoff1 <- res$log2FoldChange >= lfc & ! is.na(res$padj) & res$padj <= FDR
cutoff2 <- res$log2FoldChange <= -lfc & ! is.na(res$padj) & res$padj <= FDR
T24NPDMSO <- rownames(res[cutoffs,])
T24NPDMSO_Low <- rownames(res[cutoff2,])
T24NPDMSO_High <- rownames(res[cutoff1,])
message(sprintf("There are %s genes that are differentially expressed",sum(cutoffs)))## There are 5265 genes that are differentially expressedmessage(sprintf("There are %s genes that are induced",sum(cutoff1)))## There are 2629 genes that are inducedmessage(sprintf("There are %s genes that are repressed",sum(cutoff2)))## There are 2636 genes that are repressed6.4.2.18 MA plot
DESeq2::plotMA(res)
6.4.2.19 Volcano plot
volcanoPlot(res, lfc=1)
6.4.2.20 Heatmap
hm <- heatmap.2(t(scale(t(vsd[cutoffs,]))),col=hpal,
Colv=FALSE,dendrogram = "row",trace = "none",labRow = FALSE,
distfun = pearson.dist, labCol = paste(samples$Nutrition,samples$AZD,sep="-")[sel],
hclustfun = function(X){hclust(X,method="ward.D2")},margins = c(6,5))
6.4.2.21 T24_NP_AZD
res <- results(dds,name = "Conditions_T24_NP_AZD_vs_T0_T0_0")
data <- data.frame(rownames(res), res$log2FoldChange, res$padj)
write.table(data, file = "Conditions_T24_NP_AZD_vs_T0_T0_0.csv", sep = ";", row.names = FALSE, dec=",")
cutoffs <- abs(res$log2FoldChange) >= lfc & ! is.na(res$padj) & res$padj <= FDR
cutoff1 <- res$log2FoldChange >= lfc & ! is.na(res$padj) & res$padj <= FDR
cutoff2 <- res$log2FoldChange <= -lfc & ! is.na(res$padj) & res$padj <= FDR
T24NPAZD <- rownames(res[cutoffs,])
T24NPAZD_Low <- rownames(res[cutoff2,])
T24NPAZD_High <- rownames(res[cutoff1,])
message(sprintf("There are %s genes that are differentially expressed",sum(cutoffs)))## There are 4803 genes that are differentially expressedmessage(sprintf("There are %s genes that are induced",sum(cutoff1)))## There are 2525 genes that are inducedmessage(sprintf("There are %s genes that are repressed",sum(cutoff2)))## There are 2278 genes that are repressed6.4.2.22 MA plot
DESeq2::plotMA(res)
6.4.2.23 Volcano plot
volcanoPlot(res, lfc=1)
6.4.2.24 Heatmap
hm <- heatmap.2(t(scale(t(vsd[cutoffs,]))),col=hpal,
Colv=FALSE,dendrogram = "row",trace = "none",labRow = FALSE,
distfun = pearson.dist, labCol = paste(samples$Nutrition,samples$AZD,sep="-")[sel],
hclustfun = function(X){hclust(X,method="ward.D2")},margins = c(6,5))
7 Comparisons of GOI lists (Venn diagrams)
7.1 GOI of nutrition at T6
- All GOI
grid.draw(venn.diagram(list(T6NPSDMSO, T6NSDMSO, T6NPDMSO),
filename=NULL,
col=pal[1:3],
category.names=c("Full","No P", "No C")))
- Downregulated genes
grid.draw(venn.diagram(list(T6NPSDMSO_Low, T6NSDMSO_Low, T6NPDMSO_Low),
filename=NULL,
col=pal[1:3],
category.names=c("Full","No P", "No C")))
- Upregulated genes
grid.draw(venn.diagram(list(T6NPSDMSO_High, T6NSDMSO_High, T6NPDMSO_High),
filename=NULL,
col=pal[1:3],
category.names=c("Full","No P", "No C")))
7.2 GOI of AZD at T6
- All GOI
grid.draw(venn.diagram(list(T6NPSAZD, T6NSAZD, T6NPAZD),
filename=NULL,
col=pal[1:3],
category.names=c("Full","No P", "No C")))
- Downregulated genes
grid.draw(venn.diagram(list(T6NPSAZD_Low, T6NSAZD_Low, T6NPAZD_Low),
filename=NULL,
col=pal[1:3],
category.names=c("Full","No P", "No C")))
- Upregulated genes
grid.draw(venn.diagram(list(T6NPSAZD_High, T6NSAZD_High, T6NPAZD_High),
filename=NULL,
col=pal[1:3],
category.names=c("Full","No P", "No C")))
7.3 Nutrition and AZD interaction at T6
7.3.1 For Carbon starvation
- Downregulated genes
grid.draw(venn.diagram(list(T6NPSAZD_Low, T6NPSDMSO_Low, T6NPAZD_Low, T6NPDMSO_Low),
filename=NULL,
col=pal[1:4],
category.names=c("AZD","Ctrl","-Suc x AZD", "-Suc")))
- Upregulated genes
grid.draw(venn.diagram(list(T6NPSAZD_High, T6NPSDMSO_High, T6NPAZD_High, T6NPDMSO_High),
filename=NULL,
col=pal[1:4],
category.names=c("AZD","Ctrl","-Suc x AZD", "-Suc")))
7.3.2 For Phosphorus starvation
- Downregulated genes
grid.draw(venn.diagram(list(T6NPSAZD_Low, T6NPSDMSO_Low, T6NSAZD_Low, T6NSDMSO_Low),
filename=NULL,
col=pal[1:4],
category.names=c("AZD","Ctrl","-Phos x AZD", "-Phos")))
- Upregulated genes
grid.draw(venn.diagram(list(T6NPSAZD_High, T6NPSDMSO_High, T6NSAZD_High, T6NSDMSO_High),
filename=NULL,
col=pal[1:4],
category.names=c("AZD","Ctrl","-Phos x AZD", "-Phos")))
7.4 Nutrition and AZD interaction at T24
7.4.1 For Carbon starvation
- Downregulated genes
grid.draw(venn.diagram(list(T24NPSAZD_Low, T24NPSDMSO_Low, T24NPAZD_Low, T24NPDMSO_Low),
filename=NULL,
col=pal[1:4],
category.names=c("AZD","Ctrl","-Suc x AZD", "-Suc")))
- Upregulated genes
grid.draw(venn.diagram(list(T24NPSAZD_High, T24NPSDMSO_High, T24NPAZD_High, T24NPDMSO_High),
filename=NULL,
col=pal[1:4],
category.names=c("AZD","Ctrl","-Suc x AZD", "-Suc")))
7.4.2 For Phosphorus starvation
- Downregulated genes
grid.draw(venn.diagram(list(T24NPSAZD_Low, T24NPSDMSO_Low, T24NSAZD_Low, T24NSDMSO_Low),
filename=NULL,
col=pal[1:4],
category.names=c("AZD","Ctrl","-Phos x AZD", "-Phos")))
- Upregulated genes
grid.draw(venn.diagram(list(T24NPSAZD_High, T24NPSDMSO_High, T24NSAZD_High, T24NSDMSO_High),
filename=NULL,
col=pal[1:4],
category.names=c("AZD","Ctrl","-Phos x AZD", "-Phos")))
7.5 Comparisons of T6 of treatments to T24 of sucrose starvation
7.5.1 Comparisons to AZD treatment
- Downregulated genes
grid.draw(venn.diagram(list(T6NPSAZD_Low, T6NPSDMSO_Low, T24NPAZD_Low, T24NPDMSO_Low),
filename=NULL,
col=pal[1:4],
category.names=c("T6_AZD","T6_Ctrl","T24 -Suc x AZD", "T24 -Suc")))
- Upregulated genes
grid.draw(venn.diagram(list(T6NPSAZD_High, T6NPSDMSO_High, T24NPAZD_High, T24NPDMSO_High),
filename=NULL,
col=pal[1:4],
category.names=c("T6_AZD","T6_Ctrl","T24 -Suc x AZD", "T24 -Suc")))
7.5.2 Comparisons to phosphorus starvation
- Downregulated genes
grid.draw(venn.diagram(list(T6NSAZD_Low, T6NSDMSO_Low, T24NPAZD_Low, T24NPDMSO_Low),
filename=NULL,
col=pal[1:4],
category.names=c("T6 -Phos. x AZD","T6 -Phos.","T24 -Suc x AZD", "T24 -Suc")))
- Upregulated genes
grid.draw(venn.diagram(list(T6NSAZD_High, T6NSDMSO_High, T24NPAZD_High, T24NPDMSO_High),
filename=NULL,
col=pal[1:4],
category.names=c("T6 -Phos. x AZD","T6 -Phos.","T24 -Suc x AZD", "T24 -Suc")))
8 Export list of DEGs
write.table(T6NPSDMSO_Low, file = "T6NPSDMSO_Low.csv", sep = ";", row.names = FALSE)
write.table(T6NSDMSO_Low, file = "T6NSDMSO_Low.csv", sep = ";", row.names = FALSE)
write.table(T6NPDMSO_Low, file = "T6NPDMSO_Low.csv", sep = ";", row.names = FALSE)
write.table(T6NPSDMSO_High, file = "T6NPSDMSO_High.csv", sep = ";", row.names = FALSE)
write.table(T6NSDMSO_High, file = "T6NSDMSO_High.csv", sep = ";", row.names = FALSE)
write.table(T6NPDMSO_High, file = "T6NPDMSO_High.csv", sep = ";", row.names = FALSE)
write.table(T6NPSAZD_Low, file = "T6NPSAZD_Low.csv", sep = ";", row.names = FALSE)
write.table(T6NSAZD_Low, file = "T6NSAZD_Low.csv", sep = ";", row.names = FALSE)
write.table(T6NPAZD_Low, file = "T6NPAZD_Low.csv", sep = ";", row.names = FALSE)
write.table(T6NPSAZD_High, file = "T6NPSAZD_High.csv", sep = ";", row.names = FALSE)
write.table(T6NSAZD_High, file = "T6NSAZD_High.csv", sep = ";", row.names = FALSE)
write.table(T6NPAZD_High, file = "T6NPAZD_High.csv", sep = ";", row.names = FALSE)
write.table(T24NPSDMSO_Low, file = "T24NPSDMSO_Low.csv", sep = ";", row.names = FALSE)
write.table(T24NSDMSO_Low, file = "T24NSDMSO_Low.csv", sep = ";", row.names = FALSE)
write.table(T24NPDMSO_Low, file = "T24NPDMSO_Low.csv", sep = ";", row.names = FALSE)
write.table(T24NPSDMSO_High, file = "T24NPSDMSO_High.csv", sep = ";", row.names = FALSE)
write.table(T24NSDMSO_High, file = "T24NSDMSO_High.csv", sep = ";", row.names = FALSE)
write.table(T24NPDMSO_High, file = "T24NPDMSO_High.csv", sep = ";", row.names = FALSE)
write.table(T24NPSAZD_Low, file = "T24NPSAZD_Low.csv", sep = ";", row.names = FALSE)
write.table(T24NSAZD_Low, file = "T24NSAZD_Low.csv", sep = ";", row.names = FALSE)
write.table(T24NPAZD_Low, file = "T24NPAZD_Low.csv", sep = ";", row.names = FALSE)
write.table(T24NPSAZD_High, file = "T24NPSAZD_High.csv", sep = ";", row.names = FALSE)
write.table(T24NSAZD_High, file = "T24NSAZD_High.csv", sep = ";", row.names = FALSE)
write.table(T24NPAZD_High, file = "T24NPAZD_High.csv", sep = ";", row.names = FALSE)9 GO analysis
9.1 T6 DMSO
9.1.1 Induced genes
GO <- gopher(T6NPSDMSO_High,background=rownames(vsd)[rowSums(vsd)>0],url="athaliana",alpha=2)## Loading required package: jsonlite##
## Attaching package: 'jsonlite'## The following object is masked from 'package:purrr':
##
## flatten## No enrichments found in task: mapmanwrite.table(GO$go, file = "GO_T6NPSDMSO_High.csv", sep = ";", row.names = FALSE)
write.table(GO$kegg, file = "KEGG_T6NPSDMSO_High.csv", sep = ";", row.names = FALSE)
write.table(GO$pfam, file = "PFAM_T6NPSDMSO_High.csv", sep = ";", row.names = FALSE)
GO <- gopher(T6NSDMSO_High,background=rownames(vsd)[rowSums(vsd)>0],url="athaliana",alpha=2)## No enrichments found in task: mapmanwrite.table(GO$go, file = "GO_T6NSDMSO_High.csv", sep = ";", row.names = FALSE)
write.table(GO$kegg, file = "KEGG_T6NSDMSO_High.csv", sep = ";", row.names = FALSE)
write.table(GO$opfam, file = "PFAM_T6NSDMSO_High.csv", sep = ";", row.names = FALSE)
GO <- gopher(T6NPDMSO_High,background=rownames(vsd)[rowSums(vsd)>0],url="athaliana", alpha=2)## No enrichments found in task: mapmanwrite.table(GO$go, file = "GO_T6NPDMSO_High.csv", sep = ";", row.names = FALSE)
write.table(GO$kegg, file = "KEGG_T6NPDMSO_High.csv", sep = ";", row.names = FALSE)
write.table(GO$pfam, file = "PFAM_T6NPDMSO_High.csv", sep = ";", row.names = FALSE)9.1.2 Repressed genes
GO <- gopher(T6NPSDMSO_Low,background=rownames(vsd)[rowSums(vsd)>0],url="athaliana",alpha=2)## No enrichments found in task: mapmanwrite.table(GO$go, file = "GO_T6NPSDMSO_Low.csv", sep = ";", row.names = FALSE)
write.table(GO$kegg, file = "KEGG_T6NPSDMSO_Low.csv", sep = ";", row.names = FALSE)
write.table(GO$pfam, file = "PFAM_T6NPSDMSO_Low.csv", sep = ";", row.names = FALSE)
GO <- gopher(T6NSDMSO_Low,background=rownames(vsd)[rowSums(vsd)>0],url="athaliana", alpha=2)## No enrichments found in task: mapmanwrite.table(GO$go, file = "GO_T6NSDMSO_Low.csv", sep = ";", row.names = FALSE)
write.table(GO$kegg, file = "KEGG_T6NSDMSO_Low.csv", sep = ";", row.names = FALSE)
write.table(GO$pfam, file = "PFAM_T6NSDMSO_Low.csv", sep = ";", row.names = FALSE)
GO <- gopher(T6NPDMSO_Low,background=rownames(vsd)[rowSums(vsd)>0],url="athaliana", alpha=2)## No enrichments found in task: mapmanwrite.table(GO$go, file = "GO_T6NPDMSO_Low.csv", sep = ";", row.names = FALSE)
write.table(GO$kegg, file = "KEGG_T6NPDMSO_Low.csv", sep = ";", row.names = FALSE)
write.table(GO$pfam, file = "PFAM_T6NPDMSO_Low.csv", sep = ";", row.names = FALSE)9.2 T6 AZD
9.2.1 Induced genes
GO <- gopher(T6NPSAZD_High,background=rownames(vsd)[rowSums(vsd)>0],url="athaliana",alpha=2)## No enrichments found in task: mapmanwrite.table(GO$go, file = "GO_T6NPSAZD_High.csv", sep = ";", row.names = FALSE)
write.table(GO$kegg, file = "KEGG_T6NPSAZD_High.csv", sep = ";", row.names = FALSE)
write.table(GO$pfam, file = "PFAM_T6NPSAZD_High.csv", sep = ";", row.names = FALSE)
GO <- gopher(T6NSAZD_High,background=rownames(vsd)[rowSums(vsd)>0],url="athaliana", alpha=2)## No enrichments found in task: mapmanwrite.table(GO$go, file = "GO_T6NSAZD_High.csv", sep = ";", row.names = FALSE)
write.table(GO$kegg, file = "KEGG_T6NSAZD_High.csv", sep = ";", row.names = FALSE)
write.table(GO$pfam, file = "PFAM_T6NSAZD_High.csv", sep = ";", row.names = FALSE)
GO <- gopher(T6NPAZD_High,background=rownames(vsd)[rowSums(vsd)>0],url="athaliana", alpha=2)## No enrichments found in task: mapmanwrite.table(GO$go, file = "GO_T6NPAZD_High.csv", sep = ";", row.names = FALSE)
write.table(GO$kegg, file = "KEGG_T6NPAZD_High.csv", sep = ";", row.names = FALSE)
write.table(GO$pfam, file = "PFAM_T6NPAZD_High.csv", sep = ";", row.names = FALSE)9.2.2 Repressed genes
GO <- gopher(T6NPSAZD_Low,background=rownames(vsd)[rowSums(vsd)>0],url="athaliana", alpha=2)## No enrichments found in task: mapmanwrite.table(GO$go, file = "GO_T6NPSAZD_Low.csv", sep = ";", row.names = FALSE)
write.table(GO$kegg, file = "KEGG_T6NPSAZD_Low.csv", sep = ";", row.names = FALSE)
write.table(GO$pfam, file = "PFAM_T6NPSAZD_Low.csv", sep = ";", row.names = FALSE)
GO <- gopher(T6NSAZD_Low,background=rownames(vsd)[rowSums(vsd)>0],url="athaliana", alpha=2)## No enrichments found in task: mapmanwrite.table(GO$go, file = "GO_T6NSAZD_Low.csv", sep = ";", row.names = FALSE)
write.table(GO$pfam, file = "PFAM_T6NSAZD_Low.csv", sep = ";", row.names = FALSE)
write.table(GO$kegg, file = "KEGG_T6NSAZD_Low.csv", sep = ";", row.names = FALSE)
GO <- gopher(T6NPAZD_Low,background=rownames(vsd)[rowSums(vsd)>0],url="athaliana", alpha=2)## No enrichments found in task: mapmanwrite.table(GO$go, file = "GO_T6NPAZD_Low.csv", sep = ";", row.names = FALSE)
write.table(GO$kegg, file = "KEGG_T6NPAZD_Low.csv", sep = ";", row.names = FALSE)
write.table(GO$pfam, file = "PFAM_T6NPAZD_Low.csv", sep = ";", row.names = FALSE)9.3 T24 DMSO
9.3.1 Induced genes
GO <- gopher(T24NPSDMSO_High,background=rownames(vsd)[rowSums(vsd)>0],url="athaliana", alpha=2)## No enrichments found in task: mapmanwrite.table(GO$go, file = "GO_T24NPSDMSO_High.csv", sep = ";", row.names = FALSE)
write.table(GO$kegg, file = "KEGG_T24NPSDMSO_High.csv", sep = ";", row.names = FALSE)
write.table(GO$pfam, file = "PFAM_T24NPSDMSO_High.csv", sep = ";", row.names = FALSE)
GO <- gopher(T24NSDMSO_High,background=rownames(vsd)[rowSums(vsd)>0],url="athaliana", alpha=2)## No enrichments found in task: mapmanwrite.table(GO$go, file = "GO_T24NSDMSO_High.csv", sep = ";", row.names = FALSE)
write.table(GO$kegg, file = "KEGG_T24NSDMSO_High.csv", sep = ";", row.names = FALSE)
write.table(GO$pfam, file = "PFAM_T24NSDMSO_High.csv", sep = ";", row.names = FALSE)
GO <- gopher(T24NPDMSO_High,background=rownames(vsd)[rowSums(vsd)>0],url="athaliana", alpha=2)## No enrichments found in task: mapmanwrite.table(GO$go, file = "GO_T24NPDMSO_High.csv", sep = ";", row.names = FALSE)
write.table(GO$kegg, file = "KEGG_T24NPDMSO_High.csv", sep = ";", row.names = FALSE)
write.table(GO$pfam, file = "PFAM_T24NPDMSO_High.csv", sep = ";", row.names = FALSE)9.3.2 Repressed genes
GO <- gopher(T24NPSDMSO_Low,background=rownames(vsd)[rowSums(vsd)>0],url="athaliana", alpha=2)## No enrichments found in task: mapmanwrite.table(GO$go, file = "GO_T24NPSDMSO_Low.csv", sep = ";", row.names = FALSE)
write.table(GO$kegg, file = "KEGG_T24NPSDMSO_Low.csv", sep = ";", row.names = FALSE)
write.table(GO$pfam, file = "PFAM_T24NPSDMSO_Low.csv", sep = ";", row.names = FALSE)
GO <- gopher(T24NSDMSO_Low,background=rownames(vsd)[rowSums(vsd)>0],url="athaliana", alpha=2)## No enrichments found in task: mapmanwrite.table(GO$go, file = "GO_T24NSDMSO_Low.csv", sep = ";", row.names = FALSE)
write.table(GO$kegg, file = "KEGG_T24NSDMSO_Low.csv", sep = ";", row.names = FALSE)
write.table(GO$pfam, file = "PFAM_T24NSDMSO_Low.csv", sep = ";", row.names = FALSE)
GO <- gopher(T24NPDMSO_Low,background=rownames(vsd)[rowSums(vsd)>0],url="athaliana", alpha=2)## No enrichments found in task: mapmanwrite.table(GO$go, file = "GO_T24NPDMSO_Low.csv", sep = ";", row.names = FALSE)
write.table(GO$kegg, file = "KEGG_T24NPDMSO_Low.csv", sep = ";", row.names = FALSE)
write.table(GO$pfam, file = "PFAM_T24NPDMSO_Low.csv", sep = ";", row.names = FALSE)9.4 T24 AZD
9.4.1 Induced genes
GO <- gopher(T24NPSAZD_High,background=rownames(vsd)[rowSums(vsd)>0],url="athaliana", alpha=2)## No enrichments found in task: mapmanwrite.table(GO$go, file = "GO_T24NPSAZD_High.csv", sep = ";", row.names = FALSE)
write.table(GO$kegg, file = "KEGG_T24NPSAZD_High.csv", sep = ";", row.names = FALSE)
write.table(GO$pfam, file = "PFAM_T24NPSAZD_High.csv", sep = ";", row.names = FALSE)
GO <- gopher(T24NSAZD_High,background=rownames(vsd)[rowSums(vsd)>0],url="athaliana", alpha=2)## No enrichments found in task: mapmanwrite.table(GO$go, file = "GO_T24NSAZD_High.csv", sep = ";", row.names = FALSE)
write.table(GO$kegg, file = "KEGG_T24NSAZD_High.csv", sep = ";", row.names = FALSE)
write.table(GO$pfam, file = "PFAM_T24NSAZD_High.csv", sep = ";", row.names = FALSE)
GO <- gopher(T24NPAZD_High,background=rownames(vsd)[rowSums(vsd)>0],url="athaliana", alpha=2)## No enrichments found in task: mapmanwrite.table(GO$go, file = "GO_T24NPAZD_High.csv", sep = ";", row.names = FALSE)
write.table(GO$kegg, file = "KEGG_T24NPAZD_High.csv", sep = ";", row.names = FALSE)
write.table(GO$pfam, file = "PFAM_T24NPAZD_High.csv", sep = ";", row.names = FALSE)9.4.2 Repressed genes
GO <- gopher(T24NPSAZD_Low,background=rownames(vsd)[rowSums(vsd)>0],url="athaliana", alpha=2)## No enrichments found in task: mapmanwrite.table(GO$go, file = "GO_T24NPSAZD_Low.csv", sep = ";", row.names = FALSE)
write.table(GO$kegg, file = "KEGG_T24NPSAZD_Low.csv", sep = ";", row.names = FALSE)
write.table(GO$pfam, file = "PFAM_T24NPSAZD_Low.csv", sep = ";", row.names = FALSE)
GO <- gopher(T24NSAZD_Low,background=rownames(vsd)[rowSums(vsd)>0],url="athaliana", alpha=2)## No enrichments found in task: mapmanwrite.table(GO$go, file = "GO_T24NSAZD_Low.csv", sep = ";", row.names = FALSE)
write.table(GO$pfam, file = "PFAM_T24NSAZD_Low.csv", sep = ";", row.names = FALSE)
write.table(GO$kegg, file = "KEGG_T24NSAZD_Low.csv", sep = ";", row.names = FALSE)
GO <- gopher(T24NPAZD_Low,background=rownames(vsd)[rowSums(vsd)>0],url="athaliana", alpha=2)## No enrichments found in task: mapmanwrite.table(GO$go, file = "GO_T24NPAZD_Low.csv", sep = ";", row.names = FALSE)
write.table(GO$kegg, file = "KEGG_T24NPAZD_Low.csv", sep = ";", row.names = FALSE)
write.table(GO$pfam, file = "PFAM_T24NPAZD_Low.csv", sep = ";", row.names = FALSE)10 Session Info
## R version 3.6.1 (2019-07-05)
## Platform: x86_64-pc-linux-gnu (64-bit)
## Running under: Ubuntu 18.04.3 LTS
##
## Matrix products: default
## BLAS/LAPACK: /usr/lib/x86_64-linux-gnu/libopenblasp-r0.2.20.so
##
## locale:
## [1] LC_CTYPE=en_US.UTF-8 LC_NUMERIC=C
## [3] LC_TIME=en_US.UTF-8 LC_COLLATE=en_US.UTF-8
## [5] LC_MONETARY=en_US.UTF-8 LC_MESSAGES=en_US.UTF-8
## [7] LC_PAPER=en_US.UTF-8 LC_NAME=C
## [9] LC_ADDRESS=C LC_TELEPHONE=C
## [11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C
##
## attached base packages:
## [1] grid parallel stats4 stats graphics grDevices utils
## [8] datasets methods base
##
## other attached packages:
## [1] jsonlite_1.6 vsn_3.54.0
## [3] VennDiagram_1.6.20 futile.logger_1.4.3
## [5] tximport_1.14.0 forcats_0.4.0
## [7] stringr_1.4.0 dplyr_0.8.3
## [9] purrr_0.3.3 readr_1.3.1
## [11] tidyr_1.0.0 tibble_2.1.3
## [13] tidyverse_1.3.0 scatterplot3d_0.3-41
## [15] RColorBrewer_1.1-2 plotly_4.9.1
## [17] pander_0.6.3 magrittr_1.5
## [19] LSD_4.0-0 limma_3.42.0
## [21] hyperSpec_0.99-20180627 ggplot2_3.2.1
## [23] lattice_0.20-38 here_0.1
## [25] gplots_3.0.1.1 DESeq2_1.26.0
## [27] SummarizedExperiment_1.16.0 DelayedArray_0.12.0
## [29] BiocParallel_1.20.0 matrixStats_0.55.0
## [31] Biobase_2.46.0 GenomicRanges_1.38.0
## [33] GenomeInfoDb_1.22.0 IRanges_2.20.1
## [35] S4Vectors_0.24.0 BiocGenerics_0.32.0
## [37] data.table_1.12.6
##
## loaded via a namespace (and not attached):
## [1] colorspace_1.4-1 rprojroot_1.3-2 htmlTable_1.13.2
## [4] XVector_0.26.0 fs_1.3.1 base64enc_0.1-3
## [7] rstudioapi_0.10 affyio_1.56.0 bit64_0.9-7
## [10] AnnotationDbi_1.48.0 lubridate_1.7.4 xml2_1.2.2
## [13] splines_3.6.1 geneplotter_1.64.0 knitr_1.26
## [16] zeallot_0.1.0 Formula_1.2-3 broom_0.5.2
## [19] annotate_1.64.0 cluster_2.1.0 dbplyr_1.4.2
## [22] BiocManager_1.30.10 compiler_3.6.1 httr_1.4.1
## [25] backports_1.1.5 assertthat_0.2.1 Matrix_1.2-18
## [28] lazyeval_0.2.2 cli_1.1.0 formatR_1.7
## [31] acepack_1.4.1 htmltools_0.4.0 tools_3.6.1
## [34] affy_1.64.0 gtable_0.3.0 glue_1.3.1
## [37] GenomeInfoDbData_1.2.2 Rcpp_1.0.3 cellranger_1.1.0
## [40] vctrs_0.2.0 preprocessCore_1.48.0 gdata_2.18.0
## [43] nlme_3.1-142 xfun_0.11 rvest_0.3.5
## [46] testthat_2.3.0 lifecycle_0.1.0 gtools_3.8.1
## [49] XML_3.98-1.20 zlibbioc_1.32.0 scales_1.1.0
## [52] hms_0.5.2 lambda.r_1.2.4 curl_4.2
## [55] yaml_2.2.0 memoise_1.1.0 gridExtra_2.3
## [58] rpart_4.1-15 latticeExtra_0.6-28 stringi_1.4.3
## [61] RSQLite_2.1.2 highr_0.8 genefilter_1.68.0
## [64] checkmate_1.9.4 caTools_1.17.1.2 rlang_0.4.2
## [67] pkgconfig_2.0.3 bitops_1.0-6 evaluate_0.14
## [70] htmlwidgets_1.5.1 bit_1.1-14 tidyselect_0.2.5
## [73] R6_2.4.1 generics_0.0.2 Hmisc_4.3-0
## [76] DBI_1.0.0 pillar_1.4.2 haven_2.2.0
## [79] foreign_0.8-72 withr_2.1.2 survival_3.1-7
## [82] RCurl_1.95-4.12 nnet_7.3-12 modelr_0.1.5
## [85] crayon_1.3.4 futile.options_1.0.1 KernSmooth_2.23-16
## [88] rmarkdown_1.18 locfit_1.5-9.1 readxl_1.3.1
## [91] blob_1.2.0 reprex_0.3.0 digest_0.6.23
## [94] xtable_1.8-4 munsell_0.5.0 viridisLite_0.3.0